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Abstract. Gradual changes in exploitation, nutrient loading, etc. produce shifts 
between alternative stable states (ASS) in ecosystems which, quite often, are 
not smooth but abrupt or catastrophic. Early warnings of such catastrophic 
regime shifts are fundamental for designing management protocols for ecosystems. 
Here we study the spatial version of a popular ecological model, involving a 
logistically growing single species subject to exploitation, which is known to 
exhibit ASS. Spatial heterogeneity is introduced by a carrying capacity parameter 
varying from cell to cell in a regular lattice. Transport of biomass among cells is 
included in the form of diffusion. We investigate whether different quantities from 
statistical mechanics -like the variance, the two-point correlation function and the 
patchiness- may serve as early warnings of catastrophic phase transitions between 
the ASS. In particular, we find that the patch-size distribution follows a power 
law when the system is close to the catastrophic transition. We also provide 
links between spatial and temporal indicators and analyze how the interplay 
between diffusion and spatial heterogeneity may affect the earliness of each of 
the observables. We find that possible remedial procedures, which can be followed 
after these early signals, are more effective as the diffusion becomes lower. Finally, 
we comment on similarities and differences between these catastrophic shifts and 
paradigmatic thermodynamic phase transitions like the liquid-vapour change of 
state for a fluid like water. 
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1. Introduction 

Ecosystems are exposed to gradual change in external conditions such as climate, 
inputs of nutrients, toxic chemicals, etc. Although it is generally assumed that these 
gradual variations produce also gradual changes in the ecosystems, occasionally sudden 
catastrophic regime shifts may occur. Recent examples of ecosystems illustrating 
such changes are the shift in Caribbean coral reefs [11 [2] , shallow lakes that become 
overgrown by floating plants savannahs that are encroached suddenly by bushes 
[H ini and lakes that shift from clear to turbid [H [7] . A simple explanation for such 
drastic shifts is that the ecosystem has alternative stable states (ASS) [SJin]. In other 
words, under the same external conditions the system can be in two or more stable 
states. Hence, when subjected to a slowly changing external factor (such as climate), 
an ecosystem may show little change until it reaches a critical point where a sudden 
shift to an alternative contrasting state occurs. The presence of ASS implies that if a 
system has gone through such a state shift, it tends to remain in the new state until the 
control variable is changed back to a much lower level. This hysteresis phenomenon 
of "history dependent" alternative equilibrium states is well known in physics. 

The simplest models for describing alternative states in ecosystems correspond to 
what are known in physics parlance as mean-field (MF) models. Neglecting all spatial 
heterogeneities, these models describe the change over time of some population that 
characterizes the state of the ecosystem. These models are easy to analyze and in cases 
without significant heterogeneity their predictions are not very different from those of 
spatial models. However, in other cases the presence of a spatial dimension profoundly 
alter population dynamics or opportunities for coexistence in the real world ^lOj . In 
fact, the oversimplification of MF models casts doubt on whether the occurrence of 
an alternative stable state could be an artifact. Moreover, verifications and predictive 
power with respect to catastrophic responses to changing environmental conditions are 
still scarce for spatially extensive ecosystems. Analysis of spatially explicit models are 
relevant for other reasons. For example, to understand phenomena like clumping and 
spatial segregation in plant communities [llj . It was shown that vegetation patches, 
which have been extensively studied for arid lands [12], can be approached as a 
pattern formation phenomenon [13]- [H]. It has been hypothesized that vegetation 
patchiness could be used as a signature of imminent catastrophic shifts between 
alternative states [15] . Evidences that the patch-size distribution of vegetation follows 
a power law were later found in arid Mediterranean ecosystems |16j . This implies that 
vegetation patches were present over a wide range of size scales, thus displaying scale 
invariance. It was also found that with increasing grazing pressure, the field data 
revealed deviations from power laws. Hence, the authors proposed that this power 
law behaviour may be a warning signal for the onset of desertification. These spatial 
early warnings complement temporal ones like the variance of time series introduced 
to detect lake eutrophication [1 [TT or the impact of pollutants [T^ . 

In this work we will consider the spatial version of a general ecological model in 
terms of a logistically growing species whose consumption, loss or removal (either by 
grazing, predation or harvesting) is represented by a saturation curve [19[ 120] . The 
MF version of this model, in terms of two parameters, is known to have ASS. In 
order to take into account the spatial heterogeneity of the landscape one of the two 

II Eutrophication is an increase in nutrients leading to an enhanced growth of aquatic vegetation or 
phytoplankton and further effects including lack of oxygen and severe reductions in water quality, 
fish, and other animal populations. 
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parameters, the local parameter, is taken as dependent on the position. The other 
parameter, the global or control parameter, is taken uniform in all the system. Our 
goal is to use this framework to analyze the following questions: 

(i) How spatial heterogeneity of the environment and diffusion of matter and 
organisms affect the existence of alternative stable states. 

(ii) Whether emergent characteristic spatial patterns are really useful as early 
warnings and how they are connected with temporal signs of catastrophic shifts. 

(iii) The search for scaling laws underlaying spatial patterns and self-organization. 

We will address these issues by measuring typical observables of statistical 
mechanics, like the spatial variance, the two-point correlation function and the 
patchiness. 

This work is organized as follows. In section [2] we review the ecological MF 
model and analyze it from the point of view of Catastrophe Theory [21] . Section [3] 
is devoted to the methods used in this study and the characterization of the steady 
states reached by the system for static situations, i.e. constant values of the uniform 
parameter. In section [5] we study the dynamic case in which the uniform parameter is 
changing with time. This combination of a varying global control parameter, modelling 
a slowly changing stressor, and a local parameter, describing the heterogeneity of 
the environment, was introduced in the case of one-dimensional models in [22j . 
Besides addressing question (P, this allows to explore question ([u]), namely possible 
spatial early warnings and their connection with temporal ones. The analysis of the 
distribution of patches, the issue (|ii3), is accomplished for changing values of the control 
parameter. The usefulness of these spatial early warnings in realistic situations and 
to implement remedial actions is analyzed in section O In section [S] we compare 
how these 'flags', indicating the onset of sudden shifts, display in the ecosystem 
under consideration and in thermodynamics. The conclusions and final comments 
are presented in section [T] 

2. Mean-Field Description 

Our starting point is the population model introduced to describe grazing systems |19j 
and later used in general for several ecosystems 20_ and in particular for the case of 
the spruce budworm [231 124] . It involves a biomass density X which evolves in time 
according to: 



where r is the intrinsic per capita growth rate, K is the carrying capacity or the 
number of individuals which can be supported in a given area within natural resource 
limits, c is the maximum consumption rate and h is a half-saturation constant i.e. 
it corresponds to the value of X such that the effective consumption is half of 
the maximum consumption rate. We can rewrite H]) in terms of non-dimensional 
quantities: t' = rt, X' ~ X/h, K' = K/h and c' = c/{hr), as 



In what follows, for simplicity, we will omit the ' for the non-dimensional variables. 
The r.h.s. of ((2|) may be thought as the gradient of a potential V associated to the 




(1) 
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problem: 



V = 



dX 



K 



= -^ + ^+c(X-arctanX) (3) 



so the equilibria correspond to the roots of the first derivative of Y : 



X' 



= 



(4) 




c 

Figure 1. Folding diagrams for different values of K. 

This equation has one or three real roots (besides the trivial unstable solution 
X = 0), corresponding to one stable equilibrium state or two alternative stable 
states (separated by an unstable one). It's interesting to notice that the presence 
of alternative stable states is linked to the functional form assumed for the density 
dependent consumption. This can be modelled by different consumption functions, 
which are subdivided in three classes: linear (or Rolling type I), hyperbolic (or HoUing 
type II) and sigmoidal (or HoUing type III) [25]. Only for the sigmoidal consumption 
there occur two stable equilibria separated by an unstable one and therefore we have 
ASS. 

In figure [1] the response curve for (jl]) is depicted for different values of K. For 
K < Kc — 3"^/^ ~ 5.196 only one stable solution exists for each c. As long as 
we consider quasi stationary evolution for increasing c, the system would exhibit a 
smooth response. On the other hand, for K > Kc, the response curve is folded 
backwards at two saddle-node bifurcation points. For certain values of c the system 
can be found either in the upper or the lower stable branch. For increasing c, the 
system starts on the upper branch and varies its state smoothly until a threshold 
value is found, where a catastrophic transition to the lower branch occurs. If at this 
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point c is decreased, we would not be able to recover the state of the system before the 
transition. Instead, the system would remain on the lower branch, until we decrease 
c enough to reach another threshold value and 'jump' to the upper branch. From a 
ecological management viewpoint, it would be desirable to anticipate these transitions. 

A general formalism for treating these catastrophic regime shifts is the Elementary 
Catastrophe Theory (ECT) developed by R. Thom [21]. However, ECT works for 
static and homogeneous (MF) systems, where there is no time or spatial dependence 
of the potential. To discuss dynamics or local properties, ECT must be extended 
by incorporating some external assumptions. A change of the control parameter, 
reflecting changes of the external conditions, modifies the form of the potential. 
Therefore, as the shape of the potential changes, an original global minimum in which 
the system sits may become a metastable local minimum because other minimum 
assumes a lower value, or it even may disappear. In this case the system must jump 
from the original global minimum to the new one. ECT does not tell us when, and 
to which minimum, the jump occurs. The criterion which determines this is called a 
convention. Before discussing conventions we need to introduce two important sets of 
points in parameter space which control structural changes of the potential. 

The first of such sets of points is the bifurcation set Sb [Hj. It divides the phase 
space into two regions corresponding either to single stability or bistability of the 
system (see figure [5]) . For the (c, K) points on this curve the second derivative of the 
potential V vanishes, so the bifurcation set is given in its parametric form by: 

0=^:^^, K=^^ forx,>l (5) 

The second set of points is called the Maxwell set Sm l26] . On the Maxwell set the 
values of V at two or more stable equilibria are equal. In our case it is defined by: 



dV 
dX 



= (6) 

Xi,X2 

V{xi) = V{X2), (7) 



(see the inset of V ior K = 7.5, c = 1.91 in figure[2]). 

Sb and Sm are connected to two commonly applied criteria or conventions. 
Systems which remain in the equilibrium that they are in until it disappears are said 
to obey the delay convention. On the other hand, systems which always seek a global 
minimum of V are said to obey the Maxwell convention. Indeed these two conventions 
correspond to two extremes in a continuum of possibilities. Furthermore, real systems 
may obey either of these two conventions depending on the rate of change of the control 
parameters or on other external conditions. When the control parameters, and so the 
shape of V , change very slowly the system tends to follow the delay convention. On 
the contrary, when the control parameters change more quickly or when perturbations 
on the system are big enough, the Maxwell convention describes better the dynamics 
(more on this below). 

3. Spatial Model 

A two dimensional spatial version of the previous mean-field model is given by: 
dX{x,y;t) ( X{x,y-t) \ X{x,y-tf , 



dt V K{x,y) J '1 + X{x,y;t) 




Figure 2. Bifurcation set (solid line) with a cusp point at c = 8/3^'''^, K = Kc 
and Maxwell set (dashed). The potential V is shown for selected values of c and 
K. 



where the carrying capacity K{x,y) is a spatial heterogeneous parameter that varies 
from point to point (while the parameter c is taken as uniform) and D is the diffusion 
coefficient measuring dispersion of X in space (given in units of the intrinsic growth 
rate 1/r from section[2]). We simulated this model in a L x L regular square lattice, so 
each cell, centred at integer coordinates (i, j), can be associated with a patch of the 
ecosystem. Each cell is connected to its four nearest neighbours i. e. the von Neumann 
neighbourhood is used. To ensure numerical stability of the discretization scheme even 
for big values of the diffusion cocfhcient, the Alternating Direction Implicit Method 
|27| is used, so evolution at each time step is divided into two stages, treating implicitly 
one of the spatial coordinates at each: 



(1 + 2a)X{i,j; t +i) - aXii + t +i) - aX{i - 1, j; t +i) 



Xii,j;t) ( 1- 
a{X{i,j 



X{t,j;t) 
K{i.j) 
l-t)+X{i,j 



X{i,j;tf 



1 



X{i,j-.t)\ 
(f-2a)X(i,j;i) 



(9) 



(1 + 2a)X{i,j;t+ I) - aX{i + f, + 1) - aX{i - l,j;t+ 1) 



X(i,j;t+i) 1 



X{i,i-t+i) 



X{i,j-t- 



K{t,j) J l + X{tJ;t+^f_ 
+ a{X{i,j + l;t +1) + X{i,j - 1; t + (1 - 2a)X(z, j; t +i) 



(10) 



where a = i and d is a reduced diffusion coefficient related with D and the lattice 

o 

spacing a by = AD/a . Periodic boundary conditions (PBC) were used and L 
ranged from 100 to 800 (in fact, for different values of L in this range, no important 
differences were found). The number of time steps is typically 1000. Depending on 
the ecosystem, each time step could correspond to a day, or a month, or a year, etc. 
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The range of values for the model parameters that we use are chosen to contain 
the region of alternative stable states determined by the MF equations: the carrying 
capacity K(i, j) varies randomly from cell to cell around a fixed spatial mean (K) = 7.5 
in the interval [—Sk,Sk] where Sk = 1-0 — 2.5. Typical values for the consumption 
rate c are between 1 and 3 and for for d are between 0.1 and 5. 

3.1. Observables 

Several quantities can be measured from the time series produced by the model: 

• The spatial mean {X): 

[i and j locate each cell of the array) . 

• The spatial variance u\: 

o\ = - {Xf (12) 

• The temporal variance computed from mean values of X at different times, 
X{t), (here we take X = {X){t) ) which is defined as: 

t'=t-T \ t' = t~T J 

for temporal bins of size r (typical values for r are from 50 to 150). 

• The patchiness or cluster structure. Clusters of high (low) X are defined as 
connected regions of cells with X{i,j,t) > Xm {X{i,j,t) < Xm) where X,n is a 
threshold value. There are different criteria to define Xm, one of which is stated 
in section m 

• The two-point correlation function for pairs of cells at (ii,ji) and {12,^2), 
separated a given distance i?, which is given by: 

G2(i?) = (X(ji,ji)X(j2,j2)) - (X(zi,ji)}(X(i2, ja)) (14) 

3.2. Stable states in heterogeneous media 

In this subsection we briefiy describe the steady states reached by the system for static 
conditions, i.e. for constant values of the control parameter. The goal is to characterize 
the different alternative states, produced by eqs. Q & (fTO|) . and their corresponding 
spatial patterns and to search for scaling laws. 

In the absence of diffusion, each cell (i, j) would end up in an equilibrium value 
that is completely determined by its carrying capacity K{i,j) and the initial value of 
X at this point. So the final state of the array would be a random distribution of values 
for X (see first row of figure [3|). On the other hand, dispersion among cells allows for 
attaining global equilibrium configurations with some kind of spatial structure (second 
and third rows of figurc[3|). Notice that this structure is more noticeable as d increases. 
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c=1.75 c=1.98 




10 20 30 40 50 10 20 30 40 50 



Figure 3. A portion of 50 X 50 cells from the original 800 X 800 lattice is 
shown, grids representing the value taken by X{i,j) at each cell at equilibrium 
for (K) = 7.5, for c = 1.75 (first column) and c = 1.98 (second column). Each 
row corresponds to d = 0, d = 0.1 and d = 0.5 respectively. 



4. Alternative stable states and early warnings 

Let us now study the effect of gradually increasing stress on the system, varying c 
from 1 to 3 in steps of 2/1000. Therefore there is an important difference with the 
results presented in the previous section: now we do not let the system "thermalize" , 
i. e. each measure is performed for a different value of the control parameter c. 

We will see that some characteristics of the spatial structure may serve as early 
warnings of catastrophic shifts of the system. 

4.I. Spatial and Temporal Variance 

In figure [3] we compute {X), a\ and in terms of increasing c with {K) ~ 7.5, 
d = 0.1 and initial condition for each X{i,j) in the interval [0, (K)]. The position of 
peak for the spatial variance, ~ 2.08, is earlier than the position of peak for the 
temporal variance in nearly 110 time steps. So cr^ works better than cr^ as a warning 
signal for the upcoming transition. The reason for this is clear. When estimating the 
temporal variance one must consider past values in the time series, which correspond 
to situations where the ecosystem is far from undergoing a transition. The spatial 
variance considers only the present values, so if a signal announcing a change is present, 
it is not obscured by averaging it with data where these indications are not present. 
However, notice that when the peak in cr^ occurs, {X) has already experienced a 
decrement of almost 50% over its initial value. 

So far we have studied the shift for increasing c. Let us see what happens when 
c is decreased. In figure [5] the hysteresis cycles, yielded by these backward shifts, are 
shown for different values of d. We observe two remarkable things. First, the peak 
in cr^ is always narrower for the backward transition than in the forward transition. 
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Figure 4. {X), tr^ and ct| for d = 0.1, {K) = 7.5. The peak of cr^ occurs at 
Cm ~ 2.08 and the peak of fr^ at c ~ 2.30. 



Second, the width of the hysteresis loop decreases with d, so diffusion tends to make 
the transition more abrupt. We wiU come back to discuss the effects of diffusion in 
greater detail later in this section. 

4-2. Correlation 

The spatial variance is a particular case {R = 0) of the two-point correlation function 
(|14p . We wonder if considering R > 1 would give further information about a coming 
catastrophic shift. In figure [5] the two point correlation is depicted for R — 0,1,2,3 
(i? is measured along rows or columns of the matrix array of system's cells) . As one 
can see, the peak of the correlation for any R occurs nearly at the same value of the 
control parameter c « Cm = 2.08. 

4-3. Patchiness: Cluster structure 

In order to study the cluster structure we must define a threshold X„i as a reference 
for the grid values X{i, j). For (K) = 7.5 and d = 0.1 the maximum in cr^ is given 
at Cm — 2.08 (figure m. The value of (X) corresponding to Cm is {X)c„, — 2.89 and 
we will take it as the threshold. In the first column of figure [7] we include snapshots 
of typical patch configurations for c = Cm — 0.1, c = Cm and c = Cm + 0.1 and in the 
second column a binary representation, i.e. dark red (blue) cells correspond to cells for 
which X > {X)c^ {X < (X)c„J. The plots at the third column are the corresponding 
cluster distributions. At c = c„i the patch-size distribution follows a power law over 
two decades -with exponent 7 « —1.1 for d = 0.1 and 7 « —0.9 for d = 0.5- which 
disappears for the smaller or greater value of c. Therefore this particular distribution 
may be considered as a signature of an upcoming catastrophic shift in the system. 

4.4- The effects of diffusion 

Now we will consider the dependence on the diffusion coefficient d of the different 
introduced spatial signals. Figure [5] and figure [5] show, respectively, the variance 
and correlation for several values of d between and 5.0. Notice that the influence 
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of diffusion on and on the correlation is just the opposite. In fact for d = 
there's almost no correlation and the peak of tr^ is maximum. On the other hand, for 
d ~ 0.5 the peak in the correlation is maximum whereas the peak for is smaller 
and much more narrower. This is because we have two opposite 'forces' operating over 
the ecosystem. On the one hand its intrinsic underlying spatial heterogeneity {K = 
^ihj) ) promotes spatial fluctuations between nearest neighbours, while the diffusion 
term tends to smooth out these differences. 

The resulting spatial patterns are shown in figure [10] for Cm and different values of 
d. For low diffusion, e.g. d = 0.1, a typical configuration of the system consists in small 
patches of arbitrary different colours (that is, large global differences, measured by the 
spatial variance, and low correlation). As d increases, nearest neighbour cells group 
into larger patches or 'supercells' of the same 'color'. The result is a lower variance and 
a higher correlation. For example, the values of cr^ and G2{1) for d — 0.1 and d — 0.5 
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0^(2) 
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c 



Figure 6. Two point correlation function for different lengths, d = 0.1, (A;) = 7.5. 




10 20 30 40 50 10 20 30 40 50 10° 10' 10^ 



Figure 7. First column: A portion of 50 X 50 cells from the original 800 X 800 
lattice is shown, grids representing the value taken by X{i,j) at each cell for 
(K) = 7.5, d = 0.1. Each row corresponds to c = 1.98, c = 2.08 and c = 2.18 
respectively. Second column: same as the first for binarized data. Third column: 
Number of clusters vs. area in logarithmic scale. 
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c 



Figure 8. Spatial variance for (K) = 7.5 and different values of the discrete 
diffusion coefficient. 




"1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 

c 



Figure 9. Two point correlation function at distance R = 1 for different values 
of the discrete diffusion coefficient, (K) = 7.5. 



are: cr^ ~ 4.315 for d = 0.1 vs. 2.292 for d = 0.5 and GzCl) ^ 1-173 for d = 0.1 
vs. 1.762 for d = 0.5. Nevertheless if d increases even more, the cofor segregation is 
so strong that at this point the size of the superceUs start to decrease fowering the 
correlation. So for d = 1.0 we have: — 1.626, G2(l) — 1.458. 

5. Usefulness of the spatial early v^rarnings 

To determine the usefulness of the warning indicators presented in the previous 
section it is necessary 1) to assess their practicality and 2) if they really allow the 
implementation of corrective actions to avoid the catastrophic shift. 

5.1. Practical considerations: dealing with incomplete and noisy information 

Calculating variances over grids consisting in a large number of sites (e.g. 400x400 
or 800x800) is easy on a computer but involves a formidable task from a measuring 
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d=0.1 d=0.6 d=1.0 




10 20 30 40 50 60 70 80 10 20 30 40 50 60 70 80 lo 20 30 40 50 60 70 80 



Figure 10. A coloured grid representing the value taken by X(i,j) at each cell 
for the value of Cm = 2.08. A portion of the lattice containing 80x80 cells is 
shown. For d=0.1: a\ ~ 4.315 and G2(l) ~ 1.173. For d=0.5: aj^ ~ 2.292 and 
G2(l) ~ 1.762. For d=l.O: a\ ~ 1.626 and G2(l) ~ 1.458. 



point of view. So, in order to assess the practical difficulty of estimating cr^, we have 
performed calculations over sample grids of different sizes Lg < L. In figure fTD we 
observe that the signal does not depend qualitatively on the number of points on the 
grid that are considered to estimate aj^. In fact, even for a very small sample of 9 
points, tJ^ still exhibits a noticeable peak. Of course, the quality of the signal improves 
with the size of the sample. Additionally, since the data from real ecosystems may 
be very noisy, it is worth considering how the presence of noise alters results. So we 
assume some level of noise by adding to c a small random value belonging to some 
interval [— Jcji^c]. In figure [12] we show {X) and cr^ for Sc = 0.5. The rise of cr^ and 
the anticipation to the temporal variance are still observed. 

5.2. Possible remedial actions 

We will study the consequences of a simple remedial action consisting in immediately 
stopping the increase of the control parameter after it reaches some threshold value 
c* . In figure [13] we show the effect of keeping c constant to c* for different values of 
c* and d. For instance, if the measure is applied at the very position of the peak of 
a^, c* = c,n — 2.08 (for (K) — 7.5), its usefulness depends on the value of d. For d 
small {d — 0.1) the decay in {X{t)) stabilizes soon to a value above 2 i.e. the system 
remains in a mixed state. On the other hand, for larger values oi d {d — 0.5) the 
decay in {X(t)) continues and the ecosystem passes to the alternative state with low 
biomass, {X{t)) < 1. This figure also shows that, for d = 0.5, the remedial measure is 
effective when applied before reaches its maximum at Cm, for c* = 1.9. We checked 
that, for moderate or high diffusion (d > 0.5), this recipe of management works if c* 
is taken between the line corresponding to Sm and the right fold line of Sb (closer to 
the first than to the second one). So a possible criterion to choose c* is as the points 
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c 



Figure 11. aj^ for d = 0.1, (K) = 7.5, 5k = 2.5 calculated on lattices of size 
Ls=3 (dotted line), La=10 and La=400 (the entire lattice). 




c 



Figure 12. Same as figure|4]for c±0.5. 
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Figure 13. (X) (black) and (blue) for (K) 



7.5 in the case of a remedial 
action consisting in keeping constant the control parameter after it reaches some 
threshold value c* . The red line indicates a threshold c* coinciding with the peak 
of (T^, c* = Cm — 2.08. Full (dashed, dash-doted) curves correspond to d=0.l 
(d=0.5). The green line points a value of c* before Cm, c* = 1.9 



6. A comparison with a thermodynamic phase transition like 
liquid-vapour: from the delay to the Maxwell convention 

Catastrophes have characteristic fingerprints or 'wave flags'. Some of the standard 
catastrophe flags are: modality, sudden jumps, hysteresis and a large or anomalous 
variance [26j . These are precisely the signals we found for the considered spatial 
heterogeneous ecological model representing a species or set of species subject to 
exploitation (grazing, harvesting or predation). 

It is interesting to analyze similarities and differences with the liquid-vapour 
transition in a fluid, like water. Therefore, the biomass density X would correspond 
to the fluid density, the liquid to the high biomass density attractor and the vapour 
to the low biomass density attractor. Let us compare the above catastrophe flags for 
the fluid vs. the ecosystem: 

• Modality: The fluid is bimodal in the neighbourhood of the liquid-gas coexistence 
curve, having well deflned liquid and gas states. So this is similar in both systems. 
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Sudden Jumps: In the case of the fluid it is certainly true that sudden jumps 
occur, since there is an abrupt increase in volume when a liquid transforms into 
vapour. However, this large change in volume occurs when a slight change in the 
temperature and pressure moves the fluid from one side of the coexistence curve 
to the other. Hence, the liquid-vapour coexistence curve can be identified with 
Sm and the water changes of state obey in general the Maxwell convention. 
On the other hand, the shift in the considered model always obeys the delay 
convention: the ecosystem remains in the higher attractor (higher values of 
X) until the bifurcation set is completely traversed. However, we have seen 
in section [2] that when perturbations are big enough to allow the switching 
between equilibria on different stability branches, the systems may follow the 
Maxwell convention. Hence we will consider the effect of a sudden perturbation 
of the environment, represented here by a sharp decrease of the average carrying 
capacity (K) followed by a slow recovery. Figure [T3] illustrates this from a MF 
point of view: K is initially equal to 7.5, and for a value of the control parameter 
c = 1.68 suddenly decreases to 6. Afterwards K increases slowly in time (as 
c does) until it reaches its original value just before the system crosses Sm at 
c = 1.915. The insets show the shape of the potential V{X) just before the 
perturbation and after recovery. 




1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 2.05 2.1 2.15 2.2 



Figure 14. Variation in X produced by a global perturbation on K which 
suddenly decreases from 7.5 to 6 and slowly recovers later. The ecosystem is 
represented by a black ball before the perturbation, gray ball at intermediate step 
and white ball after recovery. Iso-K curves for K = 7.5 and 6 are depicted. The 
red arrow represents the perturbation. Insets show the shape of the potential 
V{X) just before the perturbation and after recovery. Upper right inset: path 
followed by the system under perturbation in c — K phase space. 

What happens in the case of the spatial heterogeneous and diffusive model? In 
figure[Tn]we show the evolution of the system for a completely similar perturbation 
in (K). Instead of remaining close to the initial attractor (upper branch of 
K = 7.5), the system rapidly falls to the lower branch of K = 6.0 (which 
corresponds to the minimum value of the potential V). Next it approaches more 
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slowly to the lower branch of K — 7.5 until it arrives to it for c ~ 1.915. So one 
can conclude that this type of perturbation on the system produces a change of 
convention: from delay to Maxwell. 




- ' - ' unperturbed 
perturbed 
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c 



Figure 15. The effect on (X) of a global perturbation on {K) which suddenly 
decreases from (K) = 7.5 to 6 and slowly recovers later. Thin lines represent 
iso-K curves for K = 7.5 and K = 6.0. 



• Hysteresis: In everyday situations one does not observe hysteresis in the liquid- 
gas phase transition of water: the liquid usually boils at the same temperature 
at which the vapour condenses. In other words, water changes of state obey 
in general the Maxwell convention. Nevertheless, a careful experimentalist can 
obtain an hysteresis cycle by first raising the temperature and superheating the 
liquid, and after evaporation, cooling the gas below the condensation point. 
Indeed the coexistence curve is surrounded by two spinodal lines which determine 
the limits to superheating and supersaturation. These spinodal or fold lines can 
then be identified with Sb- 

• Anomalous Variance: When a fluid condenses (boils) from its gas (liquid) to its 
liquid (gas) state, small droplets (bubbles) are formed. As a consequence, the 
variance of the volume may become large, similarly to what happens for the 
ecosystem. 

7. Conclusion 

We have analyzed a spatial ecological model whose MF version has been widely 
used to describe different relevant processes, ranging from pest outbreaks to habitat 
desertification and harvesting of aquatic plants. This model has alternative attractors, 
and is subjected to random spatial dispersion. 

For large enough values of the diffusion coefficient d, the system self-organizes 
producing characteristic spatial patterns. 

When changing the control parameter c, the transition from one attractor to 
the other is according to the delay convention. Nevertheless, is remarkable that, 
the transition occurs according to the Maxwell convention when a large enough 
perturbation is considered. This is similar to what happens in thermodynamics. In 
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general we encounter the Maxwell convention in thermodynamics: water either boils 
or condense when it reaches the saturation temperature. However, it is possible, with 
sufficient care, to superheat water or supercool steam, although any disturbance will 
produce an immediate change of phase. 

We have considered several spatial quantities both to characterize the state of 
the ecosystem and to use them as early warnings. Providing early warning signals 
is central both for management and recovery strategies of ecosystems. One of such 
observables is the spatial variance cr^ by measuring samples of X on a grid of points. 
It was found that a grid containing few points might be sufficient for the purpose of 
extracting an appropriate signal, and that a significant growth in cr^ could serve as 
an early warning of an imminent transition. This significant growth in the spatial 
variance is still observed even in the presence of moderate noise too. This is not 
surprising since noise, on the other hand, has been addressed as a promoting factor 
over persistence of alternative stable states in other ecosystems [25] . 

The spatial variance shows an advantage over the temporal one, as cr^ soars 
before than a^. The explanation for this is simple: since the former corresponds to a 
snapshot of the present state of the system while the latter includes in its computation 
data for previous times where the fiuctuations were still small. 

The origin of the rise in ax, is tied to the emergence of spatial patterns, in 
the form of patches of high/low concentration of X. We then conclude that the 
visualization of the onset of those patches, for example by aerial or satellite imaging, 
may be another indicator of the imminence of a catastrophic shift and an effective way 
of anticipating this transition. Furthermore, we found that at the very maximum of ax 
the distribution of sizes of patches becomes power law, so this particular distribution 
could serve as an early warning. Power law distribution has also been found in other 
systems as a signature of self-organization [29l |30] . 

Another observable of interest is the two-point correlation. We found that as long 
as the diffusion coefficient d increases the peak in cr^ decreases and the correlation 
increases. This dependence on d connected to the spatial patterns that emerge as d 
increases is the result of two factors that point in opposite directions: the intrinsic 
spatial heterogeneity of the ecosystem versus the dispersion or diffusion. Therefore, for 
low diffusion, a^ is the most appropriate of these two indicators to detect catastrophic 
shifts while the correlation works less well. On the other hand, for high diffusion, the 
correlation may become a more useful quantity to analyze. 

How helpful are all these warning signals in designing effective management 
protocols? Leaving aside economic considerations (which are beyond the scope of 
our analysis) this depends on different factors. For example, on the degree of diffusion 
(the size of d): The larger the diffusion the earlier the corrective action should be taken. 
For low values of d, a drastic measure, of immediately freezing the consumption rate 
c, is effective even when it has reached the value Cm at which the spatial variance is 
maximal. For larger values of d the remedial action taken at Cm can no longer avoid 
the catastrophic shift. Instead, provided there are no large perturbations, a simple 
quantitative criterion to take the remedial action is when c is over the line Sm- We 
have also seen that abrupt changes of the environmental conditions, reflected as sudden 
large variations of the parameters, can precipitate the transition to the low biomass 
catastrophic alternative state. It is worth to remark that early warning signals just 
provide a time where it is still possible to act, but not at all when the situation is 
easily reversible. 

Of course, the quantitative details of our conclusions depend on the choice of 
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parameter values employed in our model. Nevertheless, we have verified that the 
qualitative behaviour of our results do not depend strongly on those values. Rather it 
appears that our main conclusions should hold: spatial signals -variance, correlation 
and patchiness- are earlier than the temporal variance. Furthermore, spatial patterns 
formed in the process could be the fastest detectable warning that a catastrophic 
change is about to occur. Similar results have also been found for an eutrophication 
model [31] where alternative stable states are also present. 

Finally, in Ecology, as far as we know, the studies focus either on spatial early 
warnings (for instance refs. [IHIIIS]) or on temporal signals (e.g. [T71 [THl [32] ) . The 
link between these is novel. On the other hand, in the statistical physics of systems 
close to a phase transition, the connection between spatial and temporal phenomena 
-like hysteresis, critical slowing down, long range order, etc.- is well known. While the 
theory of phase transitions is well understood in thermodynamic equilibrium, its use 
in nonequilibrium systems is rather new. However, many of the fundamental concepts 
of equilibrium phase transitions - like scaling and universality- still apply in systems 
without a hermitian Hamiltonian but rather defined by transition rates, for which 
the local time-reversal symmetry is broken |33j. Moreover, in nonequilibrium systems 
a (dynamical) scaling of variables may occur even in first order transitions, when the 
order parameter jumps at the transition. This is exactly the kind of phenomenon we 
are observing for a spatial ecological model. 
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